Abstract

keywords algal-herbivore interactions, salinity, intertidal ecology

Introduction

Variations in community structure result from both the direct effects of environmental stressors on the physiology and survival of organisms and the indirect effects on the interactions between species (Underwood 1999; Dahlhoff et al. 2002; Longtin et al. 2009; Kordas et al., 2011, Diamond et al. 2017, Barton & Ives 2014). Understanding an organism’s physiological response to abiotic stress is often an important factor in explaining the distribution of single species (Miller et al. 2011; others more recent?), and traits which determine a species’ tolerance are frequently used to predict the likelihood of a species’ persistence in the face of climate change (cite). However, predictions made solely on tolerance limits to stress can often lead to misleading results, and fail to explain the distribution and abundance of species (Davis et al. 1998; Wallingford & Sorte 2019, but see Thierry et al. 2021).

Rocky intertidal shores are a unique system to study such ecological processes as they are highly dynamic in their physical environment, being at the interface of land and sea and exposed to the conditions of both with the rise and fall of the tides (Kunze et al. 2021). Fluctuations in abiotic conditions occur daily, seasonally and in the long term in response to large scale geographical and climatological processes, further influencing the dynamic nature of intertidal communities (Helmuth et al. 2002; Hsieh et al. 2005; Menge et al. 2011). In estuarine rocky intertidal ecosystems, salinity is one of the most important drivers of the performance of organisms at multiple scales of biological organisation, and thus has cascading impacts on population and community structure (Ritter et al. 2005; Schoch et al. 2006: Covernton & Harley 2020). Exposure to fresh water can induce physiological stress responses in marine gastropods, including decreased heart rate, reduced haemolymph osmolality and mortality (De Pirro et al. 1999; Chelazzi et al. 2001; Firth and Williams 2009), as well as disrupt ecological processes such as feeding, activity, reproduction and larval development rate (Cheung 1997; Zimmerman and Pechenik 1991). Similarly, decreased salinity levels have been found to reduce the survival, development and settlement of barnacle larvae, and subsequently influence adult distribution (Qiu and Qian 1999; Dineen and Hines 1994; Starczak et al. 2011). Further, hyposaline conditions inhibit the growth and photosynthetic rate of many algal species (Luo and Liu 2011; Connan and Stengel 2011; Karsten 2007), although several algal species have demonstrated a wide salinity tolerance range (Chang et al. 1999; Rath and Adhikary 2005), as well as a capacity for local adaptation to low salinities (Nygard and Ekelund 1999, 2006; Nygard and Dring 2008). As species have individualistic responses to salinity across a food web, salinity variation can have broad implications for species distribution and abundance by mediating trophic interactions (Witman and Grange 1998; Nielsen and Gosselin 2011) and can influence larger scale spatial diversity (Zacharias and Roff 2001; Hampel et al. 2009; Rubal et al. 2012).

[Paragraph: Environmental forcing is mediated by changes in the strength of interspecific interactions.]. Examples from the literature where indirect effects have a large impact on community structure. [Not sure how to transition to the next paragraph?].

Variations in estuarine salinity occur spatially, with salinity increasing with distance from the source of freshwater runoff, as well as temporally, in response to variation in river outflow (Zacharias and Roff 2001; Ysebaert and Herman 2002). Climate change is expected to impact salinity regimes in coastal regions both by intensification of the hydrological cycle that impacts precipitation patterns (IPCC year), as well as shifting the timing and strength of the Spring freshet (cite). Alterations to salinity patterns, coupled with more extreme weather events, thus has profound implications for the likelihood of future mortality events. Here, we explore the consequences of spatial variation in salinity within the Strait of Georgia on benthic algal-herbivore interactions to better understand the direct and indirect drivers of the distribution and abundance of macroalgae, molluscan grazers, and barnacles. Specifically, we consider the potential that hyposalinity has indirect effects on trophic interactions in an estuarine rocky intertidal system to drive community structure.

Because grazers, particularly limpets, are known to control high intertidal community structure in the Northeast Pacific (Cubit 1984; Farrell 1991), and because limpet species have been shown to be highly susceptible to low salinity stress compared to algae and other invertebrates (see previous examples), we predicted that observed differences in community structure between consistently high salinity environments and periodically low salinity environments would be driven by both direct effects on survivorship, and indirect trophic effects from decreased grazing. We hypothesised that limpets and other marine gastropods are disproportionately affected by the annual decrease in surface salinity experienced in the coastal areas near the mouth of the Fraser River. We further predicted that the abundance of gastropods in these periodically low salinity environments would be lower than that of the consistently high salinity areas around the Southern Gulf Islands. This reduction in gastropod grazers would lead to a reduction in grazing pressure on palatable algae, biofilms and algal spores, resulting in greater algal cover in low salinity environments. Increased algal abundance has been shown to negatively impact barnacle species (Farrell 1991), and we therefore predicted lower barnacle abundance in low salinity sites. To test these predictions, we combined laboratory salinity tolerance trials on limpets and green algae with observational surveys and manipulative limpet exclusion/inclusion experiments at three sites within each of the two salinity regions.

Materials and Methods

Study sites
The Strait of Georgia, British Columbia, presents a unique and ideal environment for studying the effects of salinity on coastal marine ecosystems. The 220 km strait is located between Vancouver Island and mainland British Columbia, and is partially isolated from the Pacific Ocean by restricted flow through narrow channels around the northern and southern tips of the island (Fig. 1). Seasonal variation in freshwater influx via the Fraser River, regularly exceeding a mean outflow rate of more than 7000 m³/s in summer months (Environment Canada 2012), causes a corresponding variation in sea surface salinity near the Fraser plume, with an annual drop from approximately 25 psu to less than 10 psu during peak discharge (Held and Harley 2009; Halverson and Pawlowicz 2011). This effect, however, declines with increasing distance from the estuary, with waters southwest of the Southern Gulf Islands maintaining salinities of 23 psu to 32 psu year round (Tully and Dodimead 1957; Halverson and Pawlowicz 2011).

All research took place on the traditional, ancestral, and unceded territory of the Musqueam, Squamish, Tsleil-Waututh and Tsawwassen nations. We conducted field studies at three sites within each of two regions with contrasting salinity regimes: West Vancouver (LS1, LS2, LS3), and the Southern Gulf Islands (HS1, HS2, & HS3), which experience consistently high salinities year-round (Fig. 1). The Southern Gulf Island (SGIs) sites are located on the southwest side of the island chain, and are not exposed to the Fraser River plume. Because of this, West Vancouver experiences reduced salinities during the summer, while the shorelines of the SGIs that do not directly border the open waters of the strait remain at consistently high salinities year-round (Fig. 2). Sea surface temperature in the two regions is comparable, ranging from 5.0 to 18.5 in West Vancouver and 6.0 to 18.5 in the SGIs (Fisheries and Oceans Canada 2009). The tidal range is greater in West Vancouver, with extreme high tides reaching 4.7 m above Canadian chart datum (approximated as the lowest astronomical tide), compared to 3.4 m in the SGIs. All sites used in this study were composed of granite rock except for HS1, which was sandstone. Areas selected for surveys and experiments were gently sloping (<40°) bedrock, and faced a variety of compass directions (see Table S1).

Transect Surveys
We conducted monthly surveys during low tide at each of the six study sites from May to August 2011. Because the tidal range differs between the two areas, we carried out surveys at the vertical height corresponding to approximately 30% submersion time. This occurs at 2.1 m in the SGIs and 3.0 m in West Vancouver. Ten metres of transect tape were laid across the rock face and eight randomly selected points were surveyed using a 25x25 cm quadrat. We counted all motile invertebrates and quantified sessile invertebrate percent cover.

Salinity Tolerance Experiments
i). Salinity tolerance and local adaptation of L. pelta
To determine whether or not the salinity tolerance of limpets is contingent on source population, we conducted the following experiment with two populations: one from a high salinity site and one from a low salinity site. We collected L. pelta, 20±5 mm in length, from HS1 in the SGIs on June 2, 2011, from a salinity of 27 psu and from LS3, West Vancouver, from a salinity of 10 psu, on June 6, 2011. Limpets from HS1 were randomly divided into eighteen 1 L Ziploc® containers with mesh walls, for a total of six limpets in each. We placed each container inside of an aquarium containing seawater at 30 psu; the salinity of the water within these aquaria was lowered by 2.5 psu per day until a salinity of 20 psu was reached. Limpets were allowed to acclimate to this salinity for ten days. We randomly divided limpets from LS3 into an additional eighteen containers, and placed all containers into aquaria containing seawater at 10 psu, increasing the salinity at increments of 2.5 psu per day to 20 psu. These limpets were allowed to acclimate to the final salinity of 20 psu for an additional six days. After the acclimation period was complete, we randomly arranged containers into eighteen aquaria, all containing seawater at 20 psu, so that each aquarium contained one container of limpets from the high salinity site and one from the low salinity site. Aquaria were randomly assigned salinity treatments of 5, 8, 11, 14, 17 and 20 psu. Aquaria were covered, provided with compressed air and placed inside of a flow-through sea water system to maintain a water temperature of 12°C. We lowered salinities at a rate of 3 psu every 30 minutes until the desired salinity was reached, and limpets remained submerged for twenty-eight days. Each day, we examined limpets for signs of mortality, including tissue damage, discolouration and rigidity; any dead limpets were removed. The experiment continued for twenty-eight days, and limpets were not fed during this time.

ii). Salinity tolerance of Lottia spp. with tidal emersion
To determine whether or not the salinity tolerance of limpets is influenced by the periodic emersion from hyposaline conditions experienced during low tides, we conducted a salinity tolerance experiment with L. pelta and L. digitalis, which incorporated a mimic of tidal exposure. See additional methods for this experiment in the Supplementary Material.

iii). Salinity Tolerance of Ulva sp.
We collected Ulva sp. from LS2 West Vancouver, from a salinity of 28 psu. Approximately 5-6 g of blot dried Ulva sp. was placed into each of sixty-four 1 L plastic bottles. Each bottle was randomly assigned a salinity treatment between 0 and 30 psu at intervals of 2.5 psu and provided with compressed air. The 0 psu treatment contained only distilled water, while all other treatments contained combinations of filtered seawater at 31 psu and dechlorinated freshwater at 0 psu. Bottles were placed inside of a flow-through sea water system to maintain a water temperature of 12°C and provided 25±5 µmol of continuous light. After three weeks, we blot dried and weighed all samples. One sample from each treatment was randomly selected to be assessed for photosynthetic efficiency using a pulse amplitude modulation (PAM) fluorometer (Jr PAM, Heinz Walz GmbH). Light intensities were altered using a 240W Fiber Optic Illuminator (6000-1, Intralux®) and screening filters. Samples were dark acclimated for one hour before quantum yields were measured by applying a saturating light pulse after reaching a steady state. We used photosynthesis vs. irradiance curves to determine maximum photosynthetic electron transport rate (ETRmax).

Herbivore Exclusion Experiment
We manually cleared all organisms at seven subsites within each of the six study sites. Each subsite included a limpet inclusion, exclusion and control plot. Inclusions and exclusions were formed by securing two copper fences, 2.5 cm high and 28.5 cm in diameter, to the rock face using Quickcrete® quick drying cement. Copper enclosures/exclosures of this type are effective barriers to limpets (Harley 2002) and partial barriers to periwinkles (Harley 2006). We marked one circular plot within each area, also 28.5 cm in diameter, with steel bolts to serve as a control. The position of control and treatments was randomised within each subsite. We did not include copper controls in this study, as previous work has shown that partial copper treatments lead to partial effects which are difficult to interpret (Johnson 1992). Ultimately, local salinity levels exceeded the lower tolerance limit of limpets in the low salinity sites, and the inclusion plots were not effective at containing L. pelta. Because of this, we only analysed exclusion and control plots. We maintained the exclusion treatments by removing limpets, along with any other grazers found inside the rings, every two weeks. We censused plots once per month during low tide, from May to August. We used a 10x10 cm quadrat to count motile invertebrates and barnacles and estimate percent cover of algae and mussels within each treatment. Salinity samples were taken at each sampling event and measured using a refractometer (S/Mill-E, Atago Inc.).

Statistical Analyses

All analyses were completed using R version 4.1.2 (R Core Team 2020).

Transect surveys
We analysed community data from the transect surveys using the vegan package, version 2.5-7. Species abundances were first relativized with a double Wisconsin transformation. This standardised species to equal maxima, then sites to equal totals, putting equal emphasis among sample units and among species. We ordinated the data with non-metric multidimensional scaling (nMDS). We then performed a permutational multivariate analysis of variance (PERMANOVA) to test the null hypothesis that the centroids of the groups are equivalent. In order to respect the dependence of sampling blocks within a site, we restricted the permutation scheme such that all quadrats along a transect were always permuted together. In doing so, 199 permutations were run on a matrix of Bray-Curtis dissimilarities. Because PERMANOVA cannot distinguish between differences in centroid location or levels of dispersion, we also used a PERMDISP test with 999 permutations, to test if the variances of the groups are different. We then conducted a Similarity of Percentages (SIMPER) analysis to investigate which species contributed most to the observed differences in salinity regions.

Lab experiments
To test whether survival of L. pelta populations originating from a low and high salinity region differ along a gradient of salinity levels, we conducted a survival analysis with the survival package, version 3.3-1 in R. We modelled the probability of survival with the Kaplan-Meier method, which is a non-parametric method to estimate survival probability from observed survival events. We then used these model fits to calculate the restricted mean survival time (RMST) for each site by population combination. We calculated a RMST for seven and 28 days.

To determine differences in net productivity of Ulva sp., we used a least-squares regression to analyse the change in biomass before and after the treatments, as well as ETRmax at the end of the experiment.

Field exclusion experiment
Due to a lack of recruitment in Spring months, and a late Summer heat wave in August 2011 that resulted in the die-off of many species, we analysed community data in the penultimate sampling point of July only. We used the same methods as described above in the transect surveys subsection. We restricted permutations similarly, keeping treatment plots within a subsite together. We then conducted a SIMPER analysis to investigate which species contribute the most to observed differences among salinity regions, grazer treatments and the interaction between the two. Prior to running the SIMPER analysis, we removed grazers from the site-species matrix, to ensure that our analysis does not identify a species that differed among treatments due to experimental manipulation. Finally, we fit generalised linear models, using the glmmTMB package - version 1.1.2.3., to the abundance of two species which were consistently identified as the most influential taxa in the SIMPER analysis, as well as had abundance patterns that were markedly altered by our experimental design: Ulva sp., and Chthamalus dalli. Both of these models included an interaction between the fixed effects, region and treatment, and site as a random effect nested within salinity region. To model percent cover of Ulva sp., we first attempted to fit a model with a beta error distribution, but this model failed to converge, therefore we used a tweedie error distribution with variance among regions modelled independently. We modelled abundance of C. dalli with a negative binomial error distribution, with a zero-inflation parameter. Both models were fitted with a type III sums of squares. We checked model diagnostics with the DHARMa package, version 0.4.5, by running K-S tests, Levene’s test and plotting scaled residuals against each predictor variable.

Results

Transect surveys
Communities belonging to low or high salinity regions had definite and substantial differences in community composition from one another. As summer progressed, communities in the low salinity region became slightly less variable. Salinity region (PERMANOVA, F = 81.18, P = 0.005; Table 1), month (PERMANOVA, pseudo-F = 4.48, p = 0.005; Table #), and the interaction between the two (PERMANOVA, pseudo-F = 1.70, p = 0.005; Table #) all had a significant impact on community composition. Dispersion among salinity groups was equal (PERMDISP, pseudo-F = 0.23, P > 0.05; table 2), although May and June have less variance than September (PERMDISP, pseudo-F = 6.52, p = 0.001; table 2). As groups with different dispersions may result in misleadingly low p values, PERMANOVA results must be interpreted cautiously. A SIMPER analysis revealed the following species contributed the most to differences between salinity groups: Mytilus trossulus, Balanus glandula, Fucus distichus, Chthamalus dalli, and the Petrocelis phase of Mastocarpus sp (table 1). Low salinity sites were composed of more M. trossulus and F. distichus, while high salinity sites were composed of more B. glandula, C. dalli, and Mastocarpus sp.

Tolerance Experiments

i). Salinity tolerance and local adaptation of L. pelta
The restricted mean survival time (RMST) for the LS3 population was higher than RMST of the HS1 population at the two lowest salinity treatments of 5 and 8 psu. Interestingly, when the RMST is calculated for the first seven days of the experiment, HS1 had a higher RMST than LS3 at a salinity of 11 psu, however this trend was reversed for an RMST calculated at 28 days (Fig. 3). Populations did not differ in their RMST in salinities above 11 psu.

ii). Salinity tolerance of Lottia spp. with tidal emersion
Survival of L. pelta was not impacted by tidal emersion, but see the Supplementary Material for more detailed results.

iii). The effect of salinity treatment on the net productivity of Ulva sp. is represented by the following second order polynomial: y = -0.915 + .342x - 0.0114x^2 (Fig. 4a, R²=0.482, P<0.001), with the greatest gain in mass at 15 psu and net losses at both 0 psu and 30 psu. ETRmax showed a similarly unimodal relationship: y = 28.2 + 4.84x - 0.13x^2 (Fig. 4b; R²=0.712, P = 0.00198), with a maximum value at 20 psu and minimum at 0 psu.

Herbivore exclusion experiment

Herbivore exclusion had little influence on community structure in the low salinity region, but had a large effect in the high salinity region. Notably, the communities in the high salinity herbivore exclusion plots were more similar to low salinity communities (with or without herbivores) than to high salinity plots with herbivores. Salinity region did not have a significant effect on community structure (PERMANOVA, pseudo-F = 0.375, P = 0.375), but treatment (PERMANOVA, pseudo-F = 3.41, p value = 0.005; Table #), as well as the interaction between salinity region and treatment (PERMANOVA, pseudo-F = 1.52, p value = 0.025; Table #) did have a significant effect. Dispersion among salinity groups was unequal (PERMDISP, F = 17.3, P = 0.001), as well as dispersion among treatment groups (PERMDISP, F = 7.5, p = 0.007).

Ulva sp., B. glandula, and C. dalli, contributed the most to differences in community composition among the high salinity with grazers and the high salinity without grazer plots (SIMPER, table #?). Similarly, the same three species also had the highest contribution to between group differences in the high salinity and low salinity plots with grazers. Excluding grazers differentially affected the abundance of C. dalli and Ulva sp. across salinity regions. Grazer excluded plots in the high salinity region had an order of magnitude more C. dalli than the control plots in the same region (grazer x salinity, x2 = 17.1, df = 1, P<0.001; figure #). Excluding grazers had the opposite impact on Ulva sp.; when grazers were absent, Ulva sp. decreased by [magnitude], but only in the high salinity sites (grazer x salinity, x2 = 14.99, df = 1, P < 0.001). There was also a consistent pattern of more Ulva sp. in the high salinity sites, ( x2 = 20.06, df = 1, P < 0.001).

Discussion

Acknowledgements

Declarations

Funding This study was funded by ?

Conflicts of interest The authors declare that they have no conflict of interest.

Ethics approval All applicable institutional and/or national guidelines for the care and use of animals were followed.

Consent to participate Not applicable

Consent for publication Not applicable

Availability of data and material The data sets generated and analysed, and the code used to analyse the data during this study are available in the [] repository.

Code availability

Authors’ contributions

CDGH initially conceived the idea for the study, CDGH, BK and TC developed the methodology, TC performed the experiments, TC and BK conducted the field work, SE and TC analyzed the data, SE and TC wrote the manuscript. All authors contributed editorial advice.

References

We would like to thank those who helped in the field; Daniel Hepler, Jocelyn Nelson, Kat Anderson, Joshua Elzam, Andres Cisneros, Kyle Demes, Rebecca Gooding, Jonathan Coyle, Steve Fuchs, Manon Picard, Gerald Singh, Jacob Uber and Darah Gibson. We would also like to thank Kyle Demes for his assistance with the PAM fluorometry experiments. This project was supported by the National Sciences and Engineering Research Council of Canada.

Tables

Figure captions

Figures

figure 1 figure 2

figure 3

figure 4

figure 5 figure 5 figure 6 figure 7